Analysis

• Run CCS with the different parameters (Testing_CCS_Parameters_Part1.sh, Testing_CCS_Parameters_Part2a.sh)
• Post CCS, run Iso-Seq3 pipeline (LIMA, REFINE, CLUSTER) to obtain the number of reads (FL, FLNC, Testing_CCS_Parameters_Part2b.sh)
• Map with ERCC, and proceed with Cupcake

Datawrangle from Cupcake output files

Use all output files from Cupcake to determine the number of FL reads associated with ERCC
1. group.txt: Output from further collapsing of polished transcripts from Iso-Seq3 to unique isoforms (denoted by PB.ID)
2. gff: Important to find which PB.ID denotes to ERCC
3. abundance: For the number of FL reads associatd with that PB.ID (note this would be the summation of FL reads across the collapsing the polished transcripts)
In theory, only need gff and abundance, however interesting to see how many FL transcripts are further collapsed

Successful CCS

Range of proporition of CCS generated can vary as much as 70% - 85%, with the main determining factor being RQ value particularly with only 1 full pass.
At 0.99 (default setting/Iso-Seq paper), number of passes do not really make a difference.

num_reads$p1

Failed CCS

The majority of polymerase reads failing to generate CCS is due to “lacking full passes” and “CCS below minimum RQ” (hence varying the two parameters).

num_reads$p2

Successful FLNC

Full-length reads are generated from CCS by removal of cDNA primers. CCS with unwanted orientations are removed and are oriented 5’ to 3’. With more stringent parameters for CCS generation, fewer reads are filtered out subsequently.

num_reads$p3

Failed FLNC

Majority of CCS reads failing to generate FL due to undesired 3p-3p pairs (cDNA synthesis?)

num_reads$p4

All Reads

Significant different in number of reads generated from RQ threshold of 0.99 and other testing RQ parameters, where number of passes does not have that much of an impact. 0.9 with 3 passes similar to 0.95 and other passes

num_reads$p5

Read Count

Similar ERCC read count observed across all parameters. However, one ERCC not detected from 0.99RQ whereas detected for other RQ. 0.99RQ is too stringent and missing real isoforms. 3 passes are clustered together for 0.8,0.9 and 0.95RQ.

draw_heatmap("Part1")

Heatmap of K18 with the FL Read count of detected ERCCs

Correlation

Strong correlation between FL count and ERCC Amount - Iso-Seq is quantitative! - No significant difference across the different parameters

# merge with ERCC conc
ERCC_output <- merge(collapsed_TOFU_reads, ERCC_conc, by.x = "V1", by.y = "ERCC_ID")

magic_for(silent = TRUE)
for (i in unique(ERCC_output$variable)) {
  df <- ERCC_output[ERCC_output$variable == i,]
  corr <- cor(df$amount_of_ERCC, df$count_fl)
  put(corr)
}
corr_output <- magic_result_as_dataframe() %>%
  mutate(num_passes = word(i, c(2), sep = fixed("_"))) %>% 
  mutate(RQ = word(i, c(4), sep = fixed("_"))) 

levels(corr_output$num_passes) <- factor(c("1","2","3"))
levels(corr_output$num_passes) <- factor(c("1","3"))
plot_line_no_reorder(corr_output, "num_passes", "corr", "Number of Passes", "Pearson's Correlation: FL Count & Amounts of ERCC","RQ") +
  mytheme + theme(legend.position = "bottom") +
  scale_colour_manual(values = c(wes_palette("Chevalier1")[3], wes_palette("Darjeeling1")[2], wes_palette("Darjeeling1")[1], wes_palette("Darjeeling2")[5])) 

Sensitivity

ERCC detected from 1 sample.
- Note, only 10% of polyermase reads were generated to make CCS
- will need to test on 1 entire sample
- does adding new sample increase ERCC detection

# which ERCC detected 
df <- ERCC_output[ERCC_output$variable == "K18_1_passes_0.8",]
for (ERCC_ID in unique(ERCC_conc$ERCC_ID)){
  Detected <- ifelse(ERCC_ID %in% noquote(df$V1), "Yes","No")
  put(Detected)
}

merge(magic_result_as_dataframe(), ERCC_conc, by = "ERCC_ID", all = T) %>% 
  ggplot(., aes(x = reorder(ERCC_ID, -amount_of_ERCC), y = log2(amount_of_ERCC),colour = Detected)) + geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  mytheme + 
  labs(x = "ERCC", y= "Log2 Stock ERCC amount")

Successful CCS

Proportion of CCS reads generated do not differ significantly between the number of passes, with sample (hence sequencing depth) being the most driving factor).

num_reads$p1

Failed CCS

The majority of polymerase reads failing to generate CCS is due to “lacking full passes” and “CCS below minimum RQ” (hence varying the two parameters).

num_reads$p2

Successful FLNC

Full-length reads are generated from CCS by removal of cDNA primers. CCS with unwanted orientations are removed and are oriented 5’ to 3’. With more stringent parameters for CCS generation, fewer reads are filtered out subsequently.

num_reads$p3

Failed FLNC

Majority of CCS reads failing to generate FL due to undesired 3p-3p pairs (cDNA synthesis?)

num_reads$p4

All Reads

The effects of using passes 1 or passes 3 remain constant in terms of the number of reads filtered through.

num_reads$p5

Read Count

Shown below is a heatmap of the FL read counts across ERCC in Samples K23 and M21 (fully generated CCS) using 0.9RQ, and minimum passes 1 and 3. Note, FL read count from each ERCC is from the summation of counts from all isoforms mapped to that ERCC (i.e with Example above PB.1.1 - PB.1.5 all mapped to ERCC-00002). Sample M21 had a significantly lower yield and number of FL reads compared to sample K23.

It appears that with great sequencing depth, the number of passes do not make a difference even in detecting lowly-expressed genes (as with Sample K23). However, with samples of lower sequencing depth (as with sample M21), a pass of 3 may generate a few false negatives.

## Heatmap of the number of FL reads associated per ERCC across the samples
# Create matrix of the heatmap for input 
reads_heatmap <- collapsed_TOFU_reads %>% filter(V3 == "transcript") %>% group_by(variable, V1) %>% tally(count_fl) %>% 
  spread(., variable, n) %>% 
  column_to_rownames(., var = "V1") %>% 
  as.matrix()

# for ERCC not detected - replace NA as 1, as required value for heatmap clustering
# note replace with 1 rather than 0 as log10 1 = 0 
reads_heatmap[is.na(reads_heatmap)] <- 1

# Create a metadata dataframe for additional information to annotate the heatmap 
reads_metadata <- as.data.frame(unique(collapsed_TOFU_reads$variable)) %>% 
  `colnames<-`(c("variable")) %>%
  mutate(num_passes = word(variable, c(2), sep = fixed("_"))) %>% 
  mutate(RQ = word(variable, c(4), sep = fixed("_"))) %>%
  mutate(Sample = word(variable, c(1), sep = fixed("_")))
# only keep the relevant columns
ann <- data.frame(reads_metadata$num_passes, reads_metadata$RQ, reads_metadata$Sample)
# determine the colours and format for the annotations
colnames(ann) <- c('Passes', 'RQ','Sample')
colours <- list('Passes' = c('1' = alpha(wes_palette("Darjeeling2")[4],0.7), '2' = wes_palette("Zissou1")[2],'3' =wes_palette("Darjeeling2")[2]),
                'RQ' = c('0.8' = wes_palette("Chevalier1")[3], '0.9' = wes_palette("Darjeeling1")[2],
                         '0.95' = wes_palette("Darjeeling1")[1],'0.99' = wes_palette("Darjeeling2")[5]), 
                'Sample' = c('K18' = wes_palette("GrandBudapest1")[2],'K23' = wes_palette("GrandBudapest1")[3], 'M21' = wes_palette("GrandBudapest1")[4]))
colAnn <- HeatmapAnnotation(df = ann, which = 'col',
                            col = colours, annotation_width = unit(c(1, 4), 'cm'), gap = unit(1, 'mm'))

# Combine matrix and annotation  
ht_list <- Heatmap(log10(reads_heatmap), 
        top_annotation = colAnn,
        show_row_dend = FALSE, show_column_names = FALSE, col=plasma(100),
        name = 'FL Count (Log10)',
        row_names_gp = gpar(fontsize = 10),
        heatmap_legend_param = list(direction = "horizontal")) 


draw(ht_list, heatmap_legend_side = "top", annotation_legend_side = "top",merge_legend = TRUE)

Correlation

High correlation (r = ~0.85) between number of FL reads and known amounts of ERCC. Of note, no difference between the different passes and the samples (despite difference in sequencing depth)

# Read in amounts of ERCC and merge with FL reads 
ERCC_conc <- read.csv("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/ONT/ERCC/ERCC_calculations.csv", header = T)[-1,] 
ERCC_output <- merge(collapsed_TOFU_reads, ERCC_conc, by.x = "V1", by.y = "ERCC_ID") %>%
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  mutate(log2_FL_reads = log2(count_fl)) 

# loop through the different datasets and aggregate the number of FL reads associated with each ERCC 
# note summation of FL reads across transcripts collapsed to ERCC
# calculate the correlation of the known amount of ERCC and summation of FL reads
magic_for(silent = TRUE)
for (i in unique(ERCC_output$variable)) {
  df <- ERCC_output[ERCC_output$variable == i & ERCC_output$V3 == "transcript",]
  df_sum <- aggregate(df$count_fl, by=list(V1=df$V1), FUN=sum) %>% left_join(ERCC_conc, by = c("V1"="ERCC_ID"))
  
  df_sum <- df_sum %>% mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>% mutate(log2_FL_reads = log2(x)) 
  corr <- cor(df_sum$log2_amount_of_ERCC, df_sum$log2_FL_reads)
  put(corr)
}

# Prepare for plots 
corr_output <- magic_result_as_dataframe() %>%
  mutate(num_passes = word(i, c(2), sep = fixed("_"))) %>% 
  mutate(RQ = word(i, c(4), sep = fixed("_"))) %>% 
  mutate(sample = word(i, c(1), sep = fixed("_")))
levels(corr_output$num_passes) <- factor(c("1","3"))

# Plot the number of passes with the correlation of FL Read count to Amount of ERCC
plot_line_no_reorder(corr_output, "num_passes", "corr", "Number of Passes", "Pearson's Correlation: Log2(FL Count) & Log2(Amounts of ERCC)","sample") +
  mytheme + theme(legend.position = "bottom") +
  scale_colour_manual(values = c(wes_palette("GrandBudapest1")[2],wes_palette("GrandBudapest1")[3],wes_palette("GrandBudapest1")[4]))

Sensitivity

Sensitivity curve of the concentration of ERCC detected per sample (with 1 pass) - no significant difference in number of ERCC detected was observed with the different passes; of note, as shown in heatmap. Unsurprisisingly, the detection of ERCC is primarily driven by sequencing depth (i.e run yield) rather than by the parameters. However, for samples with low sequencing depth, using less stringent parameters can reduce false negatives.

# sensitivity <variable>
# <variable> = name of dataset 
# Aim: draw a sensitivity curve of the number of detected ERCC by ascending known concentration 
  # 1. Extract the unique row of ERCC and known concentration (this is to remove any duplicated rows for ERCC with more than one isoform)
  # 2. Arrange ERCC concentration in ascending order
  # 3. 
sensitivity <- function(..){
  # Subset for each dataset
  df <- ERCC_output[ERCC_output$variable %in% ..,] %>% 
    # Extract the unique row of ERCC and known concentration (this is to remove any duplicated rows for ERCC with more than one isoform)
    # Arrange ERCC concentration in ascending order for cumulative observations
    select(V1,amount_of_ERCC,variable) %>% unique() %>% arrange(amount_of_ERCC) %>% 
    # Cumulative rolling count for each unique ERCC detected from lowest to highest ERCC concentration
    mutate(count = cumsum(!duplicated(V1)))
  df
}


bind_rows(sensitivity("K23_1_passes_0.9"), sensitivity("M21_1_passes_0.9"), sensitivity("K18_1_passes_0.9")) %>% 
  ggplot(.,aes(x = amount_of_ERCC, y = count, colour = variable)) + geom_point(size = (0.2)) + 
  mytheme + labs(x = "Amount of ERCC", y = "Frequency of ERCC") + theme(legend.position = c(0.8,0.2)) + geom_smooth(se = FALSE) +
  scale_colour_manual(values = c(wes_palette("GrandBudapest1")[2], wes_palette("GrandBudapest1")[3], wes_palette("GrandBudapest1")[4]),labels = c("K18","K23","M21"), name = "Sample")

The addition of samples for detecting more “novel” genes/isoforms are only important if not able to achieve initial high sequencing depth with the first sample.

bind_rows(
  sensitivity(c("K23_1_passes_0.9","M21_1_passes_0.9")) %>% select(-variable) %>% mutate(count = cumsum(!duplicated(V1))) %>% mutate(Samples = "K23+M21"), 
   sensitivity(c("K23_1_passes_0.9","K18_1_passes_0.9")) %>% select(-variable) %>% mutate(count = cumsum(!duplicated(V1))) %>% mutate(Samples = "K23+K18"),
  sensitivity("M21_1_passes_0.9") %>% select(-variable) %>% mutate(count = cumsum(!duplicated(V1))) %>% mutate(Samples = "M21"),
  sensitivity("K23_1_passes_0.9") %>% select(-variable) %>% mutate(count = cumsum(!duplicated(V1))) %>% mutate(Samples = "K23")) %>% 
  ggplot(.,aes(x = amount_of_ERCC, y = count, colour = Samples)) + geom_point(size = (0.2)) + 
  mytheme + labs(x = "Amount of ERCC", y = "Frequency of ERCC") + theme(legend.position = c(0.8,0.2)) + geom_smooth(se = FALSE) +
  scale_colour_manual(values = c(wes_palette("GrandBudapest1")[1], wes_palette("GrandBudapest1")[3], wes_palette("GrandBudapest2")[4], wes_palette("GrandBudapest1")[2]))

Group_Reads.txt

Shown below is a snapshot of group_reads.txt from TOFU cupcake. Of note, Cupcake TOFU is further collapse HQ-polished FL transcripts into unique isoforms. However, a few of these unique isoforms are mapped to the same locus region (i.e. gene with multiple isoforms), despite only one isoform known per ERCC - indicating redundancy in collapsing using default parameters from Cupcake TOFU.

head(collapsed_TOFU_grouped_reads[1:4,]) %>% kbl() %>% kable_styling() %>% 
  footnote(general = "As shown, unique isoform PB.1.4 has 2 identical transcripts (transcript/29277,transcript/29793) that are further collapsed by Cupcake TOFU. Naming system for the post-collapse isoform is PB.<loci_index>.<isoform_index>")

Input.gff

Shown below is a snapshot of input GFF

head(collapsed_TOFU_gff[collapsed_TOFU_gff$V3 == "transcript",]) %>% kbl() %>% kable_styling() %>% 
  footnote(general = "GFF is necessary to determine the corresponding PB.ID and mapped gene (in this case ERCC) from V1")

Abundance.txt

Shown below is a snapshot of abundance.txt with the summation of FL reads associated with each PB.ID - for PB1.1, 2870 full-length reads were detected across the 8 polished FL transcripts that were further collpased by Cupcake TOFU.

head(collapsed_TOFU_abundance_reads[1:4,]) %>% kbl() %>% kable_styling() %>% 
  footnote(general = "count_fl = Number of full-length reads")

Plot

# Tally the number of isoforms mapped to ERCC, and the distribution of the number of isoforms 
num_isoform <- collapsed_TOFU_reads[collapsed_TOFU_reads$V3 == "transcript",] %>% 
  group_by(V1,variable) %>% tally() %>% 
  group_by(n, variable) %>% tally() %>% 
  `colnames<-`(c("num_of_isoforms","variable","count"))

# Tally the total number of ERCCs detected across the datasets for percentage
total_num_isoform <- num_isoform %>% group_by(variable) %>% tally(count) %>%  `colnames<-`(c("variable","total_count"))
num_isoform %>%
  full_join(., total_num_isoform, by = "variable") %>% 
  mutate(perc = count/total_count * 100) %>%
  mutate(sample = word(variable, c(1), sep = fixed("_"))) %>%
  mutate(passes = word(variable, c(2), sep = fixed("_"))) %>%
  ggplot(., aes(x = num_of_isoforms, y = perc, fill = passes)) + geom_bar(stat = "identity",position = position_dodge()) +
  facet_grid(~ sample) + 
  mytheme + labs(x = "Number of Isoforms", y = "ERCC (%)") + 
  scale_fill_manual(values = c(alpha(wes_palette("Darjeeling2")[4],0.7), wes_palette("Darjeeling2")[2]))

Table

redundancy <- collapsed_TOFU_reads[collapsed_TOFU_reads$V3 == "transcript",] %>% 
  group_by(V1,variable) %>% tally() %>% select(V1,variable,n) %>% `colnames<-`(c("ERCC","Dataset","Number of isoforms"))

DT::datatable(redundancy, caption = "Number of isoforms of ERCC")

Rationale

Interestingly, all the “unique” isoforms have the same sequence, but are truncated in length. Note: the sequences highlighted are identical across all the “unique” isoforms. Shown are sequences of isoforms classified as unique by Cupcake TOFU, all mapped to ERCC-00130

Redundant Isoforms, truncated at 3’ end

Subsequently, went back to original CCS reads, and noticed all these isoforms have polyA tails. Hence, classified by Iso-Seq3 pipeline as unique isoforms after removal of polyA- also explains why not removed by Cupcake TOFU. PolyA due to random priming from oligoDT during cDNA synthesis?  the region where polyAs occur coincidentally reside with genomic As (though only 3 or 4 As)

Redundant Isoforms all with polyA

Partially-resolved redundancy

TAMA (partially) successfully removes redundant isoforms. Post-filtering with TAMA (shaded dark green), number of isoforms per ERCCs reduced significantly, with no difference between the number of passes. Note still some ERCCs with isoforms not fully removed.

TAMA_filter_dir <- "/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/WholeTranscriptome/Testing/TOFU/TAMA_FILTER"

TAMA_filtered <- list.files(path = TAMA_filter_dir, pattern = "bed", full.names = T) %>% lapply(., function(x) read.table(x) %>% select(V1,V4) %>% mutate(V1 = as.character(V1)) %>% mutate(V4 = as.character(V4)))
names(TAMA_filtered) <- list.files(path = TAMA_filter_dir, pattern = "bed")

# Tally the number of isoforms mapped to ERCC, and the distribution of the number of isoforms 
num_isoform <- TAMA_filtered %>% bind_rows(.,.id = "variable") %>%
  group_by(variable,V1) %>% tally() %>%
  group_by(n, variable) %>% tally() %>% 
  `colnames<-`(c("num_of_isoforms","variable","count"))

# Tally the total number of ERCCs detected across the datasets for percentage
total_num_isoform <- num_isoform %>% group_by(variable) %>% tally(count) %>%  `colnames<-`(c("variable","total_count"))
num_isoform %>%
  full_join(., total_num_isoform, by = "variable") %>% 
  mutate(perc = count/total_count * 100) %>%
  mutate(sample = word(variable, c(1), sep = fixed("_"))) %>%
  mutate(passes = word(variable, c(2), sep = fixed("_"))) %>%
  mutate(file = factor(word(variable, c(3), sep = fixed(".")), levels = c("bed12","bed"))) %>%
  ggplot(., aes(x = num_of_isoforms, y = perc, fill = file)) + geom_bar(stat = "identity",position = position_dodge()) +
  facet_grid(passes ~ sample) + 
  mytheme + labs(x = "Number of Isoforms", y = "ERCC (%)") + 
  scale_fill_manual(values = c(alpha(wes_palette("Darjeeling1")[2],0.3), wes_palette("Darjeeling1")[2]), 
                    labels = c("Pre-Filter","Post-Filter"), name = "") + 
  theme(legend.position = "bottom")

Retained unique ERCCs

While TAMA removes redundant isoforms, it does not remove unique ERCCs. This can be seen that the ERCCs that were detected pre-TAMA (grey-shaded) are still detected post-TAMA (blue) in that it’s not shaded in blue in the heatmap (0 isoforms), only difference is that the number of the highly expressed ERCCs isoforms are drastically reduced (> 6 isoforms to <= 2 isoforms).

## Heatmap of the number of FL reads associated per ERCC across the samples
# Create matrix of the heatmap for input 
reads_heatmap <- TAMA_filtered %>% bind_rows(.,.id = "variable") %>%  group_by(variable,V1) %>% tally() %>% 
  spread(., variable, n) %>% 
  column_to_rownames(., var = "V1") %>% 
  as.matrix()
reads_heatmap[is.na(reads_heatmap)] <- 0

# Create a metadata dataframe for additional information to annotate the heatmap 
reads_metadata <- as.data.frame(unique(total_num_isoform$variable)) %>% 
  `colnames<-`(c("variable")) %>%
  mutate(num_passes = word(variable, c(2), sep = fixed("_"))) %>% 
  mutate(RQ = word(variable, c(4), sep = fixed("_"))) %>%
  mutate(Sample = word(variable, c(1), sep = fixed("_"))) %>% 
  mutate(Filtering = factor(word(variable, c(3), sep = fixed(".")), levels = c("bed12","bed"), labels = c("Pre-Filter","Post-Filter")))
# only keep the relevant columns
ann <- data.frame(reads_metadata$num_passes, reads_metadata$Sample, reads_metadata$Filtering)
# determine the colours and format for the annotations
colnames(ann) <- c('Passes','Sample','Filtering')
colours <- list('Sample' = c('K18' = wes_palette("GrandBudapest1")[2],
                             'K23' = wes_palette("GrandBudapest1")[3], 
                             'M21' = wes_palette("GrandBudapest1")[4]),
                'Passes' = c('1' = alpha(wes_palette("Darjeeling2")[4],0.7), 
                             '2' = wes_palette("Zissou1")[2],
                             '3' =wes_palette("Darjeeling2")[2]),
                'Filtering' = c('Post-Filter' = wes_palette("Darjeeling1")[2],
                           'Pre-Filter' = wes_palette("Chevalier1")[3]))

colAnn <- HeatmapAnnotation(df = ann, which = 'col',
                            col = colours, annotation_width = unit(c(1, 4), 'cm'), gap = unit(1, 'mm'))


# Combine matrix and annotation  
ht_list <- Heatmap(reads_heatmap, 
        top_annotation = colAnn,
        show_row_dend = FALSE, show_column_names = FALSE, col=plasma(100),
        name = 'Number of Isoforms',
        row_names_gp = gpar(fontsize = 10),
        heatmap_legend_param = list(direction = "horizontal")) 


draw(ht_list, heatmap_legend_side = "top", annotation_legend_side = "top",merge_legend = TRUE)

FL read count

The FL reads per ERCC is similar pre- and post-TAMA, indicating that TAMA are removing the lowly-expressed, truncated isoforms. Note, significant difference in detecting ERCCs not determined by number of passes but sequencing depth.

## Heatmap of the number of FL reads associated per ERCC across the samples
# Create matrix of the heatmap for input 
reads_heatmap <- TAMA_filtered %>% bind_rows(.,.id = "variable") %>%
  mutate(Filtering = variable) %>%
  mutate(variable = str_remove(.$variable,".bed12")) %>% 
  mutate(variable = str_remove(.$variable,".bed")) %>% 
  mutate(PB.ID = word(.$V4, c(2), sep = fixed(";"))) %>%
  mutate(PB.ID = ifelse(!is.na(.$PB.ID), .$PB.ID, .$V4)) %>% 
  mutate(ID = paste0(variable,";", PB.ID))

reads_heatmap_final <- merge(reads_heatmap %>% select(ID, Filtering), 
                       collapsed_TOFU_reads %>% mutate(ID = paste0(variable,";", PB.ID)) %>% filter(V3 == "transcript") %>%
                         select(V1, ID, count_fl),by = "ID") %>% 
  unique(.) %>% 
  group_by(V1, Filtering) %>% tally(count_fl) %>% 
  spread(., Filtering, n) %>% 
  .[complete.cases(.[ , 1]),] %>% 
  column_to_rownames(., var = "V1") %>% 
  as.matrix()

# Create a metadata dataframe for additional information to annotate the heatmap 
reads_metadata <- as.data.frame(unique(total_num_isoform$variable)) %>% 
  `colnames<-`(c("variable")) %>%
  mutate(num_passes = word(variable, c(2), sep = fixed("_"))) %>% 
  mutate(Sample = word(variable, c(1), sep = fixed("_"))) %>% 
  mutate(Filtering = factor(word(variable, c(3), sep = fixed(".")), levels = c("bed12","bed"), labels = c("Pre-Filter","Post-Filter")))
# only keep the relevant columns
ann <- data.frame(reads_metadata$num_passes, reads_metadata$Sample, reads_metadata$Filtering)
# determine the colours and format for the annotations
colnames(ann) <- c('Passes','Sample','Filtering')
colours <- list('Sample' = c('K18' = wes_palette("GrandBudapest1")[2],
                             'K23' = wes_palette("GrandBudapest1")[3], 
                             'M21' = wes_palette("GrandBudapest1")[4]),
                'Passes' = c('1' = alpha(wes_palette("Darjeeling2")[4],0.7), 
                             '2' = wes_palette("Zissou1")[2],
                             '3' =wes_palette("Darjeeling2")[2]),
                'Filtering' = c('Post-Filter' = wes_palette("Darjeeling1")[2],
                           'Pre-Filter' = wes_palette("Chevalier1")[3]))

colAnn <- HeatmapAnnotation(df = ann, which = 'col',
                            col = colours, annotation_width = unit(c(1, 4), 'cm'), gap = unit(1, 'mm'))

  
# for ERCC not detected - replace NA as 1, as required value for heatmap clustering
# note replace with 1 rather than 0 as log10 1 = 0 
reads_heatmap_final[is.na(reads_heatmap_final)] <- 1
# Combine matrix and annotation  
ht_list <- Heatmap(log10(reads_heatmap_final), 
        top_annotation = colAnn,
        show_row_dend = FALSE, show_column_names = FALSE, col=plasma(100),
        name = 'FL Count (Log10)',
        row_names_gp = gpar(fontsize = 10),
        heatmap_legend_param = list(direction = "horizontal")) 

draw(ht_list, heatmap_legend_side = "top", annotation_legend_side = "top",merge_legend = TRUE)

SQANTI2 annotation

SQANTI2_prefilter <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/WholeTranscriptome/Testing/TOFU/TAMA_FILTER/SQANTI2/K23_3_passes_0.9_TAMAFiltered.collapsed.filtered.rep_classification.txt", header=T, as.is=T, sep="\t") %>% mutate(SQANTI2 = "Pre-Filter")
SQANTI2_postfilter <- read.table("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/WholeTranscriptome/Testing/TOFU/TAMA_FILTER/SQANTI2/K23_3_passes_0.9_TAMAFiltered.collapsed.filtered.rep_classification.filtered_lite_classification.txt", header=T, as.is=T, sep="\t") %>% mutate(SQANTI2 = "Post-Filter")

bind_rows(SQANTI2_prefilter,SQANTI2_postfilter) %>% group_by(chrom, SQANTI2) %>% tally() %>% group_by(n, SQANTI2) %>% tally() %>% ggplot(., aes(x = n, y = nn, fill = SQANTI2)) + geom_bar(stat = "identity", position = position_dodge()) + 
  mytheme + labs(x = "Number of ERCCs", y = "Number of isoforms")

SQANTI2_post_filter_counts <- SQANTI2_postfilter %>% group_by(chrom) %>% tally(FL)
ERCC_conc <- read.csv("/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/ONT/ERCC/ERCC_calculations.csv", header = T)[-1,] 
ERCC_output <- merge(SQANTI2_post_filter_counts, ERCC_conc, by.x = "chrom", by.y = "ERCC_ID") %>% 
  mutate(log2_amount_of_ERCC = log2(amount_of_ERCC)) %>%
  mutate(log2_FL_reads = log2(n)) 

density_plot(ERCC_output,"log2_amount_of_ERCC","log2_FL_reads", "Amount of ERCC (Log2)", "Number of FL reads (Log2)","")

## 
##  Pearson's product-moment correlation
## 
## data:  dat[[x.var]] and dat[[y.var]]
## t = 22.451, df = 42, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9288241 0.9785393
## sample estimates:
##      cor 
## 0.960772 
## 
## [1] ""
## [1] "Correlation betweenlog2_amount_of_ERCCandlog2_FL_reads"
## [1] "corr.value0.960771963049452"
## [1] "p.value5.13793254372933e-25"

Partially-resolved redundancy

Re-running Cupcake significantly reduced redundacy, though note still some ERCCs with 2 isoforms. No difference in terms of number of isoforms collapsed between passes.

# Input output files generated from re-running Cupcake  
working_tofu_dir <- "/gpfs/mrc0/projects/Research_Project-MRC148213/sl693/WholeTranscriptome/Testing/TOFU/RETOFU"
TOFU_read("collapsed.group.txt", "Testing_grouped_reads")
TOFU_read("collapsed.gff", "Testing_TOFU_gff")
TOFU_read("abundance.txt", "Testing_TOFU_abundance_reads")

# merge for full information 
Testing_TOFU_reads <- Testing_grouped_reads %>% 
  full_join(., Testing_TOFU_gff, by = c("PB.ID","variable")) %>% 
  full_join(., Testing_TOFU_abundance_reads, by = c("PB.ID","variable")) %>% 
  mutate(threediff = word(variable, c(5), sep = fixed("_")))

# Tally the number of isoforms mapped to ERCC, and the distribution of the number of isoforms 
num_isoform <- Testing_TOFU_reads[Testing_TOFU_reads$V3 == "transcript",] %>% 
  group_by(V1, variable) %>% tally() %>% 
  group_by(n, variable) %>% tally()  %>% 
  `colnames<-`(c("num_of_isoforms","variable","count"))

# Tally the total number of ERCCs detected across the datasets for percentage
total_num_isoform <- num_isoform %>% group_by(variable) %>% tally(count) %>%  `colnames<-`(c("variable","total_count"))
num_isoform %>%
  full_join(., total_num_isoform, by = "variable") %>% 
  mutate(perc = count/total_count * 100) %>%
  mutate(sample = word(variable, c(1), sep = fixed("_"))) %>%
  mutate(passes = word(variable, c(2), sep = fixed("_"))) %>%
  mutate(num_of_isoforms = as.factor(num_of_isoforms)) %>%
  ggplot(., aes(x = num_of_isoforms, y = perc, fill = passes)) + geom_bar(stat = "identity",position = position_dodge()) +
  facet_grid(~ sample) + 
  mytheme + labs(x = "Number of Isoforms", y = "ERCC (%)") + theme(legend.position = "bottom") +
  scale_fill_manual(values = c(alpha(wes_palette("Darjeeling2")[4],0.7), wes_palette("Darjeeling2")[2]))

FL read count

To re-collapse with Cupcake, the abundance file is generated collapsing the output fasta from the first run with Cupcake, such as the collapse is with “PB.19.1|ERCC-00096:1-1091(+)|transcript/30027” rather than “transcript/24540,transcript/25156…”. The abundance file therefore does not accurately reflect that from the clustering from Iso-Seq3 pipeline.

Given that it’s quite circular, and the abundance output is not appropriate, this solution is not taken forward.

## Heatmap of the number of FL reads associated per ERCC across the samples
# Create matrix of the heatmap for input 

First_Second_TOFU <- bind_rows(collapsed_TOFU_reads %>% mutate(variable = paste0(variable,"_First")), 
                           Testing_TOFU_reads %>% mutate(variable = paste0(variable,"_Second")))

reads_heatmap <- First_Second_TOFU %>% filter(V3 == "transcript") %>% group_by(variable, V1) %>% tally(count_fl) %>% 
  spread(., variable, n) %>% 
  column_to_rownames(., var = "V1") %>% 
  as.matrix()

# for ERCC not detected - replace NA as 1, as required value for heatmap clustering
# note replace with 1 rather than 0 as log10 1 = 0 
reads_heatmap[is.na(reads_heatmap)] <- 1

# Combine matrix and annotation  
Heatmap(log10(reads_heatmap), show_row_dend = FALSE, col=plasma(400), name = 'FL Count (Log10)')